Show code
# Create two scenarios with different steepness values
scenario_low <- create_scenario(h = 0.6, "Low,h=0.6")
scenario_high <- create_scenario(h = 0.9, "High, h=0.9")
scenarios <- list(scenario_low, scenario_high)Bering Sea Pollock-like Population Simulation test
Comparing Low and High Steepness Scenarios with constant HCR
This document is a draft for discussion purposes only.
1 Introduction
Recruitment productivity under climate change is uncertain, and stock-recruitment steepness (\(h\)) is a key parameter that drives sustainable harvest levels. Evidence from climate-ecosystem modeling and recruitment studies suggests steepness and productivity can shift with warming and prey dynamics (Hollowed et al. 2020; Holsman et al. 2020; Spencer et al. 2019; Szuwalski et al. 2023). Because steepness informs reference points and F proxies, mis-specification can change expected yields and risk profiles (Szuwalski and Punt 2025; Punt et al. 2024).
At the same time, there is pressure to incorporate environmental covariates into harvest control rules (HCRs) to improve adaptability. These approaches can help track productivity changes, but they also add interpretive uncertainty and require strong validation (Punt et al. 2024; Szuwalski et al. 2023). In contrast, tier-based HCRs remain transparent and easier to communicate, with clear links from biomass to fishing mortality (North Pacific Fishery Management Council 2024). This study asks whether a simple, transparent HCR can be robust to plausible productivity shifts without relying on covariates.
We use a stylized pollock-like, age-structured model and compare two steepness scenarios (h=0.6 and h=0.9) (North Pacific Fishery Management Council 2024). The HCR is held constant (F40% SPR with a biomass ramp from 0.4\(B_0\) to 0.05\(B_0\)), consistent with the Tier 3 framework used for Bering Sea groundfish (North Pacific Fishery Management Council 2024; Ianelli et al. 2011). The goal is to evaluate whether a single rule performs reasonably across both low and high productivity regimes, and to quantify the tradeoffs if productivity assumptions are wrong (Szuwalski and Punt 2025).
This two-scenario comparison also frames a plausible climate-driven trajectory: productivity could shift over time from a high-steepness regime to a low-steepness regime (or vice versa) as environmental conditions change. Interpreting the scenarios as endpoints of a time-varying process clarifies the management problem: the HCR is calibrated under one regime but must remain acceptable if the system transitions to the other. The results therefore inform both static mis-specification risk and the dynamic case of a climate-driven regime shift.
Key design choices include: - 15 age classes with Bering Sea-like growth and maturity schedules (Ianelli et al. 2011) - Beverton-Holt recruitment with lognormal variability (CV 0.7) (Spencer et al. 2019) - 50-year projections with 100 Monte Carlo replicates per scenario - No explicit environmental covariates; steepness shifts serve as a reduced-form proxy for climate-driven productivity changes (Hollowed et al. 2020; Holsman et al. 2020)
Outputs include reference points, time series of biomass and catch, and sensitivity comparisons across steepness values. We move from model setup to reference points in Table 1, then visualize recruitment and yield tradeoffs in Figure 1 and Figure 2. Simulation outcomes are summarized in Figure 4 and Figure 5, with sensitivity results in Figure 7, similar to HCR robustness evaluations in recent studies (Punt et al. 2024).
2 Model Description
This section documents the equations and inputs used to generate the reference points and simulations shown later.
2.1 Population Dynamics
We define the annual accounting that drives the time series trajectories shown in Figure 4.
2.1.1 Numbers at Age
The population is tracked as numbers-at-age \(N_{a,t}\) where \(a\) denotes age and \(t\) denotes year. The population evolves according to:
Recruitment (age 1): \[N_{1,t} = R_t\]
Ages 2 through 14: \[N_{a,t} = N_{a-1,t-1} \cdot e^{-Z_{a-1,t-1}}\]
Plus group (age 15): \[N_{15,t} = N_{14,t-1} \cdot e^{-Z_{14,t-1}} + N_{15,t-1} \cdot e^{-Z_{15,t-1}}\]
where \(Z_{a,t} = M + F_t \cdot s_a\) is the total instantaneous mortality rate.
These equations determine the age structure that feeds into biomass and catch trajectories in Figure 4.
2.2 Biological Parameters
Biological inputs control growth, maturity, and vulnerability, which in turn shape the yield curves in Figure 2 and the SSB trends in Figure 4.
2.2.1 Weight-at-Age
Weight-at-age is specified empirically based on the updated schedule:
| Age | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15+ |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Weight (kg) | 0.029 | 0.180 | 0.372 | 0.491 | 0.613 | 0.749 | 0.883 | 1.014 | 1.133 | 1.243 | 1.348 | 1.388 | 1.469 | 1.553 | 1.791 |
These weights translate numbers-at-age into biomass and catch, which are summarized in Figure 4 and Figure 2.
2.2.2 Maturity and Selectivity
Maturity is specified directly using fixed proportions at age, while fishery selectivity follows a logistic function:
| Age | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15+ |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Proportion mature | 0.000 | 0.008 | 0.289 | 0.641 | 0.842 | 0.901 | 0.947 | 0.963 | 0.970 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 |
\[s_a = \frac{1}{1 + e^{-\gamma_s(a - a_{50}^s)}} \quad \text{with } a_{50}^s = 4.0\]
Maturity governs SSB while selectivity determines fishing impact, which together shape the trajectories in Figure 4.
2.2.3 Spawning Stock Biomass
SSB is defined as female spawning biomass only:
\[SSB_t = \sum_{a=1}^{15} N_{a,t} \cdot W_a \cdot m_a \cdot 0.5\]
This definition is used directly in the SSB time series shown in Figure 4 and the Monte Carlo envelope in Figure 5.
2.3 Stock-Recruitment Relationship
Recruitment dynamics are summarized visually in Figure 1 and influence the variability seen in Figure 5.
2.3.1 Beverton-Holt Model with Steepness
The expected recruitment follows the Beverton-Holt stock-recruitment relationship:
\[R_t = \frac{4hR_0 \cdot SSB_{t-1}}{SSB_0(1-h) + SSB_{t-1}(5h-1)}\]
Where steepness (\(h\)) represents the fraction of \(R_0\) produced when SSB is at 20% of unfished levels. Higher steepness indicates a more resilient stock that maintains recruitment even at lower biomass levels.
As shown in Figure 1, this formulation differs across the two steepness scenarios.
2.3.2 Recruitment Variability
Recruitment includes multiplicative log-normal process error with CV = 0.7:
\[R_t = \bar{R}_t \cdot e^{\epsilon_t - \sigma_R^2/2}\]
where \(\epsilon_t \sim N(0, \sigma_R^2)\) and \(\sigma_R = \sqrt{\ln(1 + 0.7^2)} \approx 0.614\).
This stochasticity drives the spread in Figure 5 and contributes to the variability summarized in Table 4.
3 Scenario Setup
We instantiate the two steepness scenarios used throughout the tables and figures.
These objects feed the reference points in Table 1 and the graphical comparisons in Figure 1 and Figure 2.
All simulations initialize the population at the scenario’s Fmsy equilibrium rather than the unfished B0 state.
4 Reference Point Comparison
We summarize equilibrium reference points for each scenario to quantify differences in Fmsy, Bmsy, and MSY (see Table 1).
Show code
comparison_df <- compare_reference_points(scenarios)
knitr::kable(comparison_df, align = c("l", "c", "c", "c", "c", "c", "c", "c"))| Scenario | Steepness | Fmsy | Bmsy (Mt) | MSY (Mt) | SSB0 (Mt) | Bmsy/SSB0 | R_msy (B) |
|---|---|---|---|---|---|---|---|
| Low,h=0.6 | 0.6 | 0.3517 | 2.66 | 1.31 | 8.29 | 0.321 | 9.98 |
| High, h=0.9 | 0.9 | 1.2805 | 1.59 | 2.00 | 8.29 | 0.192 | 12.08 |
Table 1 provides the numeric basis for the contrasts discussed below and for the vertical markers in Figure 2.
4.1 Key Differences
Show code
ref_low <- scenario_low$ref_points
ref_high <- scenario_high$ref_pointsWe expand on the table summary by translating the reference point differences into percent changes.
The contrast between scenarios reveals important management implications:
| Metric | Low Steepness (h=0.6) | High Steepness (h=0.9) | Change |
|---|---|---|---|
| \(F_{MSY}\) | 0.352 | 1.28 | 264% higher |
| \(B_{MSY}\) (Mt) | 2.66 | 1.59 | -40% lower |
| MSY (Mt) | 1.31 | 2 | 53% higher |
| \(B_{MSY}/SSB_0\) | 0.321 | 0.192 | -40% lower |
Interpretation:
- Higher steepness allows higher fishing mortality: With \(h = 0.9\), the stock can sustain \(F_{MSY}\) = 1.28 compared to only 0.35 for \(h = 0.6\).
- Higher steepness produces higher MSY: The more productive stock yields 2 million tonnes vs. 1.31 million tonnes.
- Higher steepness requires lower \(B_{MSY}\): The stock can be fished down further while still achieving MSY.
5 Graphical Comparisons
These figures visualize how the two steepness settings diverge in recruitment, yield, and harvest control (see Figure 1, Figure 2, and Figure 3).
5.1 Stock-Recruitment Curves
As shown in Figure 1, the stock-recruitment curves differ and highlight the different Bmsy locations.
Show code
plot_sr_comparison(scenarios, colors = c("blue", "red"))
The stock-recruitment curves illustrate how steepness affects the relationship between spawning biomass and recruitment:
- Low steepness (blue): Recruitment declines more steeply as SSB decreases, making the stock more vulnerable to overfishing.
- High steepness (red): Recruitment remains relatively high even at reduced SSB, indicating greater resilience.
5.2 Yield Curves
As shown in Figure 2, equilibrium yield responds differently to fishing mortality across scenarios.
Show code
plot_yield_comparison(scenarios, colors = c("blue", "red"))
The yield curves show:
- Low steepness (blue): Yield peaks at a lower \(F\) and declines more rapidly with overfishing.
- High steepness (red): Yield peaks at a much higher \(F\) and the curve is flatter, providing more buffer against recruitment overfishing.
5.3 Harvest Control Rule (F40% with biomass ramp)
We add a sloping harvest control rule based on \(F_{40\%}\) (SPR) and unfished biomass \(B_0\). The target F values are reported in Table 2 and the rule shape is shown in Figure 3.
\[ F_t = \begin{cases} F_{40\%} & \text{if } SSB_t \ge 0.4 B_0 \\ F_{40\%} \cdot \dfrac{SSB_t - 0.05 B_0}{0.4 B_0 - 0.05 B_0} & \text{if } 0.05 B_0 < SSB_t < 0.4 B_0 \\ 0 & \text{if } SSB_t \le 0.05 B_0 \end{cases} \]
Show code
F40_low <- calc_F_spr(scenario_low$params, spr_target = 0.4)
F40_high <- calc_F_spr(scenario_high$params, spr_target = 0.4)
f40_table <- data.frame(
Scenario = c("Low steepness (h=0.6)", "High steepness (h=0.9)"),
`F40%` = round(c(F40_low, F40_high), 3),
check.names = FALSE
)
knitr::kable(f40_table)| Scenario | F40% |
|---|---|
| Low steepness (h=0.6) | 0.412 |
| High steepness (h=0.9) | 0.412 |
Show code
ssb_ratio <- seq(0, 1, by = 0.01)
F_hcr_low <- calc_sloping_hcr(ssb_ratio * ref_low$SSB0, ref_low$SSB0, F40_low)
F_hcr_high <- calc_sloping_hcr(ssb_ratio * ref_high$SSB0, ref_high$SSB0, F40_high)
y_max <- max(c(F_hcr_low, F_hcr_high)) * 1.1
hcr_df <- data.frame(
SSB_ratio = rep(ssb_ratio, 2),
F = c(F_hcr_low, F_hcr_high),
Scenario = rep(c("h=0.6 (F40%)", "h=0.9 (F40%)"), each = length(ssb_ratio))
)
hcr_colors <- c("h=0.6 (F40%)" = "blue", "h=0.9 (F40%)" = "red")
p_hcr <- ggplot(hcr_df, aes(x = SSB_ratio, y = F, color = Scenario)) +
geom_line(size = 1) +
geom_vline(xintercept = c(0.05, 0.4), linetype = "dotted", color = "gray50") +
coord_cartesian(ylim = c(0, y_max)) +
labs(x = "SSB / B0", y = "Fishing Mortality (F)", title = "Harvest Control Rule") +
scale_color_manual(values = hcr_colors) +
ggthemes::theme_few()
p_hcr
6 Simulation Results
We first compare single realizations at Fmsy (Figure 4), then summarize means (Table 3) and variability across Monte Carlo runs (Figure 5 and Table 4).
6.1 Single Simulation Comparison
As shown in Figure 4, we plot a single 50-year trajectory for each scenario at its Fmsy.
Show code
# Run simulations
sim_low <- run_simulation(scenario_low$params, F_series = ref_low$Fmsy, seed = 42)
sim_high <- run_simulation(scenario_high$params, F_series = ref_high$Fmsy, seed = 42)
years <- 1:50
make_long <- function(sim, scenario_label) {
data.frame(
Year = rep(years, 4),
Scenario = scenario_label,
Metric = rep(c("SSB (Mt)", "Recruitment (B)", "Catch (Mt)", "F"), each = length(years)),
Value = c(sim$SSB / 1e9, sim$Recruitment / 1e9, sim$Catch / 1e9, sim$F)
)
}
sim_long <- rbind(
make_long(sim_low, "h=0.6"),
make_long(sim_high, "h=0.9")
)
sim_long$Metric <- factor(
sim_long$Metric,
levels = c("SSB (Mt)", "Recruitment (B)", "Catch (Mt)", "F")
)
ref_lines <- data.frame(
Metric = c("SSB (Mt)", "SSB (Mt)", "Catch (Mt)", "Catch (Mt)"),
Scenario = c("h=0.6", "h=0.9", "h=0.6", "h=0.9"),
Value = c(ref_low$Bmsy / 1e9, ref_high$Bmsy / 1e9, ref_low$MSY / 1e9, ref_high$MSY / 1e9)
)
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
p_single <- ggplot(sim_long, aes(x = Year, y = Value, color = Scenario)) +
geom_line(size = 1) +
geom_hline(
data = ref_lines,
aes(yintercept = Value, color = Scenario),
linetype = "dashed",
size = 0.8,
inherit.aes = FALSE,
show.legend = FALSE
) +
facet_wrap(~ Metric, scales = "free_y", ncol = 2) +
scale_color_manual(values = scenario_colors) +
labs(x = "Year", y = NULL, title = "50-year simulations at Fmsy") +
ggthemes::theme_few()
p_single
6.2 Simulation Summary Statistics
Table 3 aggregates the single-simulation outcomes into mean and CV metrics.
Show code
summary_df <- data.frame(
Scenario = c("Low, h=0.6", "High, h=0.9"),
`Mean SSB (Mt)` = c(round(mean(sim_low$SSB)/1e9, 2), round(mean(sim_high$SSB)/1e9, 2)),
`CV SSB` = c(round(sd(sim_low$SSB)/mean(sim_low$SSB), 3),
round(sd(sim_high$SSB)/mean(sim_high$SSB), 3)),
`Mean Catch (Mt)` = c(round(mean(sim_low$Catch)/1e9, 2), round(mean(sim_high$Catch)/1e9, 2)),
`Mean Rec (B)` = c(round(mean(sim_low$Recruitment)/1e9, 2),
round(mean(sim_high$Recruitment)/1e9, 2)),
check.names = FALSE
)
knitr::kable(summary_df)| Scenario | Mean SSB (Mt) | CV SSB | Mean Catch (Mt) | Mean Rec (B) |
|---|---|---|---|---|
| Low, h=0.6 | 2.75 | 0.352 | 1.35 | 10.09 |
| High, h=0.9 | 1.60 | 0.413 | 2.01 | 12.13 |
6.3 Monte Carlo Simulations
As shown in Figure 5, stochastic variability across 100 runs is summarized with metrics reported in Table 4.
Show code
# Run Monte Carlo simulations
n_sims <- 100
n_years <- 50
ssb_low <- matrix(NA, n_years, n_sims)
ssb_high <- matrix(NA, n_years, n_sims)
catch_low <- matrix(NA, n_years, n_sims)
catch_high <- matrix(NA, n_years, n_sims)
for (i in 1:n_sims) {
sim_l <- run_simulation(scenario_low$params, F_series = ref_low$Fmsy, seed = i)
sim_h <- run_simulation(scenario_high$params, F_series = ref_high$Fmsy, seed = i)
ssb_low[, i] <- sim_l$SSB
ssb_high[, i] <- sim_h$SSB
catch_low[, i] <- sim_l$Catch
catch_high[, i] <- sim_h$Catch
}
ssb_low_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
SSB = as.vector(ssb_low) / 1e9,
Scenario = "h=0.6"
)
ssb_high_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
SSB = as.vector(ssb_high) / 1e9,
Scenario = "h=0.9"
)
ssb_df <- rbind(ssb_low_df, ssb_high_df)
median_df <- rbind(
data.frame(Year = 1:n_years, SSB = apply(ssb_low, 1, median) / 1e9, Scenario = "h=0.6"),
data.frame(Year = 1:n_years, SSB = apply(ssb_high, 1, median) / 1e9, Scenario = "h=0.9")
)
bmsy_df <- data.frame(
Scenario = c("h=0.6", "h=0.9"),
Bmsy = c(ref_low$Bmsy, ref_high$Bmsy) / 1e9
)
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
p_mc <- ggplot(ssb_df, aes(x = Year, y = SSB, group = interaction(Scenario, Sim), color = Scenario)) +
geom_line(alpha = 0.1, size = 0.4, show.legend = FALSE) +
geom_line(
data = median_df,
aes(x = Year, y = SSB, color = Scenario),
size = 1.1,
inherit.aes = FALSE
) +
geom_hline(
data = bmsy_df,
aes(yintercept = Bmsy, color = Scenario),
linetype = "dashed",
size = 0.8,
inherit.aes = FALSE
) +
facet_wrap(~ Scenario, ncol = 2, scales = "free_y") +
scale_color_manual(values = scenario_colors) +
labs(x = "Year", y = "Female SSB (Mt)", title = "Monte Carlo SSB trajectories") +
ggthemes::theme_few() +
guides(color = "none")
p_mc
Show code
catch_low_hcr <- matrix(NA, n_years, n_sims)
catch_high_hcr <- matrix(NA, n_years, n_sims)
for (i in 1:n_sims) {
sim_l_hcr <- run_simulation_hcr(
scenario_low$params,
F_target = F40_low,
B0 = ref_low$SSB0,
catch_cap = NULL,
seed = 2000 + i
)
sim_h_hcr <- run_simulation_hcr(
scenario_high$params,
F_target = F40_high,
B0 = ref_high$SSB0,
catch_cap = NULL,
seed = 3000 + i
)
catch_low_hcr[, i] <- sim_l_hcr$Catch
catch_high_hcr[, i] <- sim_h_hcr$Catch
}
catch_low_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
Catch = as.vector(catch_low) / 1e9,
Scenario = "h=0.6",
Policy = "Const. F"
)
catch_high_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
Catch = as.vector(catch_high) / 1e9,
Scenario = "h=0.9",
Policy = "Const. F"
)
catch_low_hcr_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
Catch = as.vector(catch_low_hcr) / 1e9,
Scenario = "h=0.6",
Policy = "FSPR HCR"
)
catch_high_hcr_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
Catch = as.vector(catch_high_hcr) / 1e9,
Scenario = "h=0.9",
Policy = "FSPR HCR"
)
catch_df <- rbind(catch_low_df, catch_high_df, catch_low_hcr_df, catch_high_hcr_df)
median_catch_df <- rbind(
data.frame(
Year = 1:n_years,
Catch = apply(catch_low, 1, median) / 1e9,
Scenario = "h=0.6",
Policy = "Const. F"
),
data.frame(
Year = 1:n_years,
Catch = apply(catch_high, 1, median) / 1e9,
Scenario = "h=0.9",
Policy = "Const. F"
),
data.frame(
Year = 1:n_years,
Catch = apply(catch_low_hcr, 1, median) / 1e9,
Scenario = "h=0.6",
Policy = "FSPR HCR"
),
data.frame(
Year = 1:n_years,
Catch = apply(catch_high_hcr, 1, median) / 1e9,
Scenario = "h=0.9",
Policy = "FSPR HCR"
)
)
policy_colors <- c("Const. F" = "#1F77B4", "FSPR HCR" = "#D95F02")
p_catch <- ggplot(
catch_df,
aes(x = Year, y = Catch, group = interaction(Scenario, Policy, Sim), color = Policy)
) +
geom_line(alpha = 0.08, size = 0.4, show.legend = FALSE) +
geom_line(
data = median_catch_df,
aes(x = Year, y = Catch, color = Policy),
size = 1.1,
inherit.aes = FALSE
) +
facet_wrap(~ Scenario, ncol = 2) +
scale_color_manual(values = policy_colors) +
labs(x = "Year", y = "Catch (Mt)", title = "Monte Carlo catch trajectories") +
ggthemes::theme_few()
p_catch
Show code
mc_stats <- data.frame(
Scenario = c("Low, h=0.6", "High, h=0.9"),
`Mean terminal SSB (Mt)` = c(round(mean(ssb_low[n_years,])/1e9, 2),
round(mean(ssb_high[n_years,])/1e9, 2)),
`P(SSB < Bmsy) terminal` = c(round(mean(ssb_low[n_years,] < ref_low$Bmsy), 2),
round(mean(ssb_high[n_years,] < ref_high$Bmsy), 2)),
`Mean catch (Mt)` = c(round(mean(catch_low)/1e9, 2),
round(mean(catch_high)/1e9, 2)),
`CV mean catch` = c(round(sd(colMeans(catch_low))/mean(catch_low), 3),
round(sd(colMeans(catch_high))/mean(catch_high), 3)),
check.names = FALSE
)
knitr::kable(mc_stats)| Scenario | Mean terminal SSB (Mt) | P(SSB < Bmsy) terminal | Mean catch (Mt) | CV mean catch |
|---|---|---|---|---|
| Low, h=0.6 | 2.58 | 0.56 | 1.26 | 0.134 |
| High, h=0.9 | 1.64 | 0.54 | 1.96 | 0.105 |
7 Fishing Mortality Sensitivity
To stress-test harvest rates, we run each scenario at constant fishing mortality values spanning 0.5-1.5 x \(F_{MSY}\). Table 5 lists the constant-F results alongside the FSPR HCR outcomes (with and without the 1.45 Mt cap), and Figure 7 visualizes the tradeoffs.
Show code
F_multipliers <- c(0.5, 0.75, 1.0, 1.25, 1.5)
n_sims_F <- 100
f_low <- run_F_sensitivity(scenario_low, F_multipliers, n_sims = n_sims_F, seed = 200)
f_high <- run_F_sensitivity(scenario_high, F_multipliers, n_sims = n_sims_F, seed = 500)
f_summary <- rbind(f_low, f_high)
f_summary$Policy <- "Const. F"
run_hcr_summary_table <- function(scenario, F_target, B0, n_sims, seed_base, catch_cap, policy_label) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
ssb_mean <- numeric(n_sims)
ssb_terminal <- numeric(n_sims)
catch_mean <- numeric(n_sims)
f_mean <- numeric(n_sims)
ssb_lt_bmsy <- logical(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation_hcr(
params,
F_target = F_target,
B0 = B0,
catch_cap = catch_cap,
seed = seed_base + i
)
ssb_mean[i] <- mean(sim$SSB)
ssb_terminal[i] <- sim$SSB[n_years]
catch_mean[i] <- mean(sim$Catch)
f_mean[i] <- mean(sim$F)
ssb_lt_bmsy[i] <- sim$SSB[n_years] < ref_points$Bmsy
}
data.frame(
Scenario = scenario$name,
Steepness = params$h,
Policy = policy_label,
F_multiplier = mean(f_mean) / ref_points$Fmsy,
F_applied = mean(f_mean),
Mean_SSB = mean(ssb_mean),
Terminal_SSB = mean(ssb_terminal),
P_SSB_lt_Bmsy = mean(ssb_lt_bmsy),
Mean_Catch = mean(catch_mean)
)
}
catch_cap <- 1.45e9
hcr_low <- run_hcr_summary_table(
scenario_low,
F40_low,
ref_low$SSB0,
n_sims_F,
seed_base = 800,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_high <- run_hcr_summary_table(
scenario_high,
F40_high,
ref_high$SSB0,
n_sims_F,
seed_base = 900,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_low_cap <- run_hcr_summary_table(
scenario_low,
F40_low,
ref_low$SSB0,
n_sims_F,
seed_base = 1000,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_high_cap <- run_hcr_summary_table(
scenario_high,
F40_high,
ref_high$SSB0,
n_sims_F,
seed_base = 1100,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_summary <- rbind(hcr_low, hcr_high, hcr_low_cap, hcr_high_cap)
f_summary_all <- rbind(f_summary, hcr_summary)
f_table <- data.frame(
Scenario = f_summary_all$Scenario,
Policy = f_summary_all$Policy,
`F multiplier` = signif(f_summary_all$F_multiplier, 2),
`F applied` = round(f_summary_all$F_applied, 3),
`Mean SSB (Mt)` = round(f_summary_all$Mean_SSB / 1e9, 2),
`Terminal SSB (Mt)` = round(f_summary_all$Terminal_SSB / 1e9, 2),
`P(SSB < Bmsy) terminal` = round(f_summary_all$P_SSB_lt_Bmsy, 2),
`Mean Catch (Mt)` = round(f_summary_all$Mean_Catch / 1e9, 2),
check.names = FALSE
)
knitr::kable(f_table)| Scenario | Policy | F multiplier | F applied | Mean SSB (Mt) | Terminal SSB (Mt) | P(SSB < Bmsy) terminal | Mean Catch (Mt) |
|---|---|---|---|---|---|---|---|
| Low,h=0.6 | Const. F | 0.50 | 0.176 | 4.11 | 4.32 | 0.08 | 1.10 |
| Low,h=0.6 | Const. F | 0.75 | 0.264 | 3.22 | 3.30 | 0.29 | 1.24 |
| Low,h=0.6 | Const. F | 1.00 | 0.352 | 2.67 | 2.50 | 0.61 | 1.31 |
| Low,h=0.6 | Const. F | 1.20 | 0.440 | 2.19 | 2.17 | 0.76 | 1.29 |
| Low,h=0.6 | Const. F | 1.50 | 0.528 | 1.87 | 1.72 | 0.90 | 1.27 |
| High, h=0.9 | Const. F | 0.50 | 0.640 | 2.37 | 2.52 | 0.08 | 1.85 |
| High, h=0.9 | Const. F | 0.75 | 0.960 | 1.85 | 1.90 | 0.45 | 1.93 |
| High, h=0.9 | Const. F | 1.00 | 1.280 | 1.55 | 1.45 | 0.67 | 1.95 |
| High, h=0.9 | Const. F | 1.20 | 1.601 | 1.39 | 1.29 | 0.77 | 2.01 |
| High, h=0.9 | Const. F | 1.50 | 1.921 | 1.21 | 1.19 | 0.89 | 1.96 |
| Low,h=0.6 | FSPR HCR | 0.90 | 0.318 | 2.72 | 2.69 | 0.57 | 1.26 |
| High, h=0.9 | FSPR HCR | 0.28 | 0.354 | 3.22 | 3.25 | 0.00 | 1.63 |
| Low,h=0.6 | FSPR HCR cap | 0.78 | 0.275 | 3.13 | 3.20 | 0.35 | 1.21 |
| High, h=0.9 | FSPR HCR cap | 0.19 | 0.237 | 4.22 | 4.40 | 0.00 | 1.35 |
Show code
plot_data <- f_summary
plot_data$Policy <- "Const. F"
plot_data$Scenario_plot <- ifelse(
plot_data$Scenario == scenario_low$name,
"h=0.6",
"h=0.9"
)
catch_cap <- 1.45e9
run_hcr_summary <- function(scenario, F_target, B0, n_sims, seed_base, catch_cap, policy_label) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
ssb_mean <- numeric(n_sims)
ssb_terminal <- numeric(n_sims)
catch_mean <- numeric(n_sims)
f_mean <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation_hcr(
params,
F_target = F_target,
B0 = B0,
catch_cap = catch_cap,
seed = seed_base + i
)
ssb_mean[i] <- mean(sim$SSB)
ssb_terminal[i] <- sim$SSB[n_years]
catch_mean[i] <- mean(sim$Catch)
f_mean[i] <- mean(sim$F)
}
data.frame(
Scenario = scenario$name,
F_multiplier = mean(f_mean) / ref_points$Fmsy,
Mean_SSB = mean(ssb_mean),
Terminal_SSB = mean(ssb_terminal),
Mean_Catch = mean(catch_mean),
Policy = policy_label
)
}
hcr_low <- run_hcr_summary(
scenario_low,
F40_low,
ref_low$SSB0,
n_sims_F,
seed_base = 800,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_high <- run_hcr_summary(
scenario_high,
F40_high,
ref_high$SSB0,
n_sims_F,
seed_base = 900,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_low_cap <- run_hcr_summary(
scenario_low,
F40_low,
ref_low$SSB0,
n_sims_F,
seed_base = 1000,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_high_cap <- run_hcr_summary(
scenario_high,
F40_high,
ref_high$SSB0,
n_sims_F,
seed_base = 1100,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_summary <- rbind(hcr_low, hcr_high, hcr_low_cap, hcr_high_cap)
hcr_summary$Scenario_plot <- ifelse(
hcr_summary$Scenario == scenario_low$name,
"h=0.6",
"h=0.9"
)
plot_long <- rbind(
data.frame(
F_multiplier = plot_data$F_multiplier,
Scenario = plot_data$Scenario_plot,
Metric = "Mean Catch (Mt)",
Value = plot_data$Mean_Catch / 1e9,
Policy = plot_data$Policy
),
data.frame(
F_multiplier = plot_data$F_multiplier,
Scenario = plot_data$Scenario_plot,
Metric = "Terminal SSB (Mt)",
Value = plot_data$Terminal_SSB / 1e9,
Policy = plot_data$Policy
),
data.frame(
F_multiplier = hcr_summary$F_multiplier,
Scenario = hcr_summary$Scenario_plot,
Metric = "Mean Catch (Mt)",
Value = hcr_summary$Mean_Catch / 1e9,
Policy = hcr_summary$Policy
),
data.frame(
F_multiplier = hcr_summary$F_multiplier,
Scenario = hcr_summary$Scenario_plot,
Metric = "Terminal SSB (Mt)",
Value = hcr_summary$Terminal_SSB / 1e9,
Policy = hcr_summary$Policy
)
)
plot_long$Metric <- factor(plot_long$Metric, levels = c("Mean Catch (Mt)", "Terminal SSB (Mt)"))
bmsy_df <- data.frame(
Scenario = c("h=0.6", "h=0.9"),
Metric = "Terminal SSB (Mt)",
Value = c(ref_low$Bmsy, ref_high$Bmsy) / 1e9
)
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
p_f <- ggplot() +
geom_line(
data = plot_long[plot_long$Policy == "Const. F", ],
aes(x = F_multiplier, y = Value, color = Scenario),
size = 1
) +
geom_point(
data = plot_long[plot_long$Policy == "Const. F", ],
aes(x = F_multiplier, y = Value, color = Scenario),
size = 2
) +
geom_point(
data = plot_long[plot_long$Policy != "Const. F", ],
aes(x = F_multiplier, y = Value, color = Scenario, shape = Policy),
size = 3,
stroke = 0.8
) +
geom_hline(
data = bmsy_df,
aes(yintercept = Value, color = Scenario),
linetype = "dashed",
size = 0.8,
inherit.aes = FALSE,
show.legend = FALSE
) +
facet_wrap(~ Metric, scales = "free_y", ncol = 2) +
scale_color_manual(values = scenario_colors) +
scale_shape_manual(values = c("FSPR HCR" = 17, "FSPR HCR cap" = 15)) +
labs(x = "F multiplier (relative to Fmsy)", y = NULL, title = "F sensitivity tradeoffs") +
ggthemes::theme_few()
p_f
Show code
n_sims_points <- n_sims_F
catch_cap <- 1.45e9
sim_constant_points <- function(scenario, f_multiplier, seed_base, n_sims) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
mean_catch <- numeric(n_sims)
terminal_ssb <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation(params, F_series = ref_points$Fmsy * f_multiplier, seed = seed_base + i)
mean_catch[i] <- mean(sim$Catch)
terminal_ssb[i] <- sim$SSB[n_years]
}
data.frame(
Scenario = scenario$name,
F_multiplier = rep(f_multiplier, n_sims),
Mean_Catch = mean_catch,
Terminal_SSB = terminal_ssb,
Policy = "Const. F"
)
}
sim_hcr_points <- function(scenario, F_target, B0, seed_base, n_sims, catch_cap, policy_label) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
mean_catch <- numeric(n_sims)
terminal_ssb <- numeric(n_sims)
mean_f <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation_hcr(
params,
F_target = F_target,
B0 = B0,
catch_cap = catch_cap,
seed = seed_base + i
)
mean_catch[i] <- mean(sim$Catch)
terminal_ssb[i] <- sim$SSB[n_years]
mean_f[i] <- mean(sim$F)
}
data.frame(
Scenario = scenario$name,
F_multiplier = mean_f / ref_points$Fmsy,
Mean_Catch = mean_catch,
Terminal_SSB = terminal_ssb,
Policy = policy_label
)
}
const_low_pts <- sim_constant_points(scenario_low, 1.0, seed_base = 1200, n_sims_points)
const_high_pts <- sim_constant_points(scenario_high, 1.0, seed_base = 1300, n_sims_points)
hcr_low_pts <- sim_hcr_points(
scenario_low,
F40_low,
ref_low$SSB0,
seed_base = 1400,
n_sims_points,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_high_pts <- sim_hcr_points(
scenario_high,
F40_high,
ref_high$SSB0,
seed_base = 1500,
n_sims_points,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_low_cap_pts <- sim_hcr_points(
scenario_low,
F40_low,
ref_low$SSB0,
seed_base = 1600,
n_sims_points,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_high_cap_pts <- sim_hcr_points(
scenario_high,
F40_high,
ref_high$SSB0,
seed_base = 1700,
n_sims_points,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
points_df <- rbind(const_low_pts, const_high_pts, hcr_low_pts, hcr_high_pts, hcr_low_cap_pts, hcr_high_cap_pts)
points_df$Scenario_plot <- ifelse(
points_df$Scenario == scenario_low$name,
"h=0.6",
"h=0.9"
)
bmsy_lookup <- ifelse(
points_df$Scenario_plot == "h=0.6",
ref_low$Bmsy,
ref_high$Bmsy
)
plot_points <- data.frame(
Catch_Mt = points_df$Mean_Catch / 1e9,
Terminal_SSB_Bmsy = points_df$Terminal_SSB / bmsy_lookup,
Scenario = points_df$Scenario_plot,
Policy = points_df$Policy,
PolicyGroup = ifelse(points_df$Policy == "Const. F", "Const. F", "SPR-based HCR")
)
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
p_fp <- ggplot(plot_points, aes(x = Catch_Mt, y = Terminal_SSB_Bmsy, color = Scenario, shape = Policy)) +
geom_point(alpha = 0.6, size = 2) +
geom_smooth(
data = plot_points,
aes(x = Catch_Mt, y = Terminal_SSB_Bmsy, linetype = PolicyGroup, group = PolicyGroup),
method = "loess",
se = FALSE,
size = 0.8,
color = "black",
inherit.aes = FALSE
) +
expand_limits(y = 0) +
scale_color_manual(values = scenario_colors) +
scale_linetype_manual(values = c("Const. F" = "solid", "SPR-based HCR" = "dashed")) +
scale_shape_manual(values = c("Const. F" = 16, "FSPR HCR" = 17, "FSPR HCR cap" = 15)) +
labs(x = "Mean Catch (Mt)", y = "Terminal SSB / Bmsy", title = "F sensitivity: individual simulation draws") +
ggthemes::theme_few()
p_fp
Show code
catch_cap <- 1.45e9
last10_stats_constant <- function(scenario, f_multiplier, seed_base, n_sims) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
last_years <- (n_years - 9):n_years
mean_catch <- numeric(n_sims)
cv_catch <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation(params, F_series = ref_points$Fmsy * f_multiplier, seed = seed_base + i)
catch_last <- sim$Catch[last_years]
mean_catch[i] <- mean(catch_last)
cv_catch[i] <- sd(catch_last) / mean(catch_last)
}
data.frame(
Scenario = scenario$name,
Mean_Catch = mean_catch,
Catch_CV = cv_catch,
Policy = "Const. F"
)
}
last10_stats_hcr <- function(scenario, F_target, B0, seed_base, n_sims, catch_cap, policy_label) {
params <- scenario$params
n_years <- params$n_years
last_years <- (n_years - 9):n_years
mean_catch <- numeric(n_sims)
cv_catch <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation_hcr(
params,
F_target = F_target,
B0 = B0,
catch_cap = catch_cap,
seed = seed_base + i
)
catch_last <- sim$Catch[last_years]
mean_catch[i] <- mean(catch_last)
cv_catch[i] <- sd(catch_last) / mean(catch_last)
}
data.frame(
Scenario = scenario$name,
Mean_Catch = mean_catch,
Catch_CV = cv_catch,
Policy = policy_label
)
}
const_low_var <- last10_stats_constant(scenario_low, 1.0, seed_base = 1600, n_sims_points)
const_high_var <- last10_stats_constant(scenario_high, 1.0, seed_base = 1700, n_sims_points)
hcr_low_var <- last10_stats_hcr(
scenario_low,
F40_low,
ref_low$SSB0,
seed_base = 1800,
n_sims_points,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_high_var <- last10_stats_hcr(
scenario_high,
F40_high,
ref_high$SSB0,
seed_base = 1900,
n_sims_points,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_low_cap_var <- last10_stats_hcr(
scenario_low,
F40_low,
ref_low$SSB0,
seed_base = 2000,
n_sims_points,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_high_cap_var <- last10_stats_hcr(
scenario_high,
F40_high,
ref_high$SSB0,
seed_base = 2100,
n_sims_points,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
var_df <- rbind(const_low_var, const_high_var, hcr_low_var, hcr_high_var, hcr_low_cap_var, hcr_high_cap_var)
var_df$Scenario_plot <- ifelse(
var_df$Scenario == scenario_low$name,
"h=0.6",
"h=0.9"
)
var_df$Mean_Catch_Mt <- var_df$Mean_Catch / 1e9
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
var_summary <- aggregate(
cbind(Mean_Catch_Mt, Catch_CV) ~ Scenario_plot + Policy,
data = var_df,
FUN = mean
)
p_var <- ggplot(var_summary, aes(x = Mean_Catch_Mt, y = Catch_CV, color = Scenario_plot, shape = Policy)) +
geom_point(size = 5) +
expand_limits(y = 0) +
scale_color_manual(values = scenario_colors) +
scale_shape_manual(values = c("Const. F" = 16, "FSPR HCR" = 17, "FSPR HCR cap" = 15)) +
labs(x = "Mean Catch (Mt, last 10 years)", y = "Catch CV (last 10 years)", title = "Catch variability vs mean catch") +
ggthemes::theme_few()
p_var
8 Alternating Steepness Simulation
8.1 Base-recruitment variability
To represent a climate-driven regime shift, we run a separate simulation where steepness alternates between high and low values every 10 years. The capped SPR HCR (F40% with a 1.45 Mt catch cap) is applied throughout, and we summarize productivity using log(recruitment/SSB) for each simulation run.
Show code
n_years_alt <- scenario_low$params$n_years
block_years <- 10
block_id <- rep(seq_len(ceiling(n_years_alt / block_years)), each = block_years)[1:n_years_alt]
h_series <- ifelse(block_id %% 2 == 1, scenario_high$params$h, scenario_low$params$h)
n_sims_alt <- 100
catch_cap_alt <- 1.45e9
F_target_alt <- calc_F_spr(scenario_low$params, spr_target = 0.4)
B0_alt <- calc_SSB0(scenario_low$params)
log_rps_list <- vector("list", n_sims_alt)
for (i in 1:n_sims_alt) {
sim_alt <- run_simulation_hcr_variable_h(
scenario_low$params,
h_series = h_series,
F_target = F_target_alt,
B0 = B0_alt,
catch_cap = catch_cap_alt,
seed = 4000 + i
)
rps <- sim_alt$Recruitment[-1] / pmax(sim_alt$SSB[-length(sim_alt$SSB)], 1e-12)
log_rps_list[[i]] <- data.frame(
Year = 2:n_years_alt,
Sim = i,
Log_R_SSB = log(rps)
)
}
log_rps_df <- do.call(rbind, log_rps_list)
median_log_rps <- aggregate(Log_R_SSB ~ Year, data = log_rps_df, FUN = median)
p_alt <- ggplot(log_rps_df, aes(x = Year, y = Log_R_SSB, group = Sim)) +
geom_line(alpha = 0.2, size = 0.4, color = "steelblue") +
geom_line(
data = median_log_rps,
aes(x = Year, y = Log_R_SSB),
inherit.aes = FALSE,
color = "black",
size = 1.1
) +
labs(x = "Year", y = "log(Recruitment / SSB)", title = "Alternating steepness (10-year blocks)") +
ggthemes::theme_few()
p_alt
8.2 Alternating Steepness with Low Recruitment Variability
We repeat the alternating-steepness experiment with low recruitment variability (\(\sigma_R = 0.1\)) to show how reduced process error changes the log(recruitment/SSB) dynamics under the capped HCR.
Show code
n_years_alt_low <- scenario_low$params$n_years
block_years_low <- 10
block_id_low <- rep(seq_len(ceiling(n_years_alt_low / block_years_low)), each = block_years_low)[1:n_years_alt_low]
h_series_low <- ifelse(block_id_low %% 2 == 1, scenario_high$params$h, scenario_low$params$h)
params_low_sigma <- scenario_low$params
params_low_sigma$rec_sigma <- 0.1
params_low_sigma$rec_cv <- sqrt(exp(params_low_sigma$rec_sigma^2) - 1)
n_sims_alt_low <- 100
catch_cap_alt_low <- 1.45e9
F_target_alt_low <- calc_F_spr(params_low_sigma, spr_target = 0.4)
B0_alt_low <- calc_SSB0(params_low_sigma)
log_rps_list_low <- vector("list", n_sims_alt_low)
for (i in 1:n_sims_alt_low) {
sim_alt_low <- run_simulation_hcr_variable_h(
params_low_sigma,
h_series = h_series_low,
F_target = F_target_alt_low,
B0 = B0_alt_low,
catch_cap = catch_cap_alt_low,
seed = 5000 + i
)
rps_low <- sim_alt_low$Recruitment[-1] / pmax(sim_alt_low$SSB[-length(sim_alt_low$SSB)], 1e-12)
log_rps_list_low[[i]] <- data.frame(
Year = 2:n_years_alt_low,
Sim = i,
Log_R_SSB = log(rps_low)
)
}
log_rps_low_df <- do.call(rbind, log_rps_list_low)
median_log_rps_low <- aggregate(Log_R_SSB ~ Year, data = log_rps_low_df, FUN = median)
p_alt_low <- ggplot(log_rps_low_df, aes(x = Year, y = Log_R_SSB, group = Sim)) +
geom_line(alpha = 0.2, size = 0.4, color = "steelblue") +
geom_line(
data = median_log_rps_low,
aes(x = Year, y = Log_R_SSB),
inherit.aes = FALSE,
color = "black",
size = 1.1
) +
labs(
x = "Year",
y = "log(Recruitment / SSB)",
title = "Alternating steepness (10-year blocks), sigma_R = 0.1"
) +
ggthemes::theme_few()
p_alt_low
8.2.1 Alternating Steepness with Regime-Specific Fmsy
As a contrast to the capped HCR, we apply the “correct” \(F_{MSY}\) for each steepness block (high or low) while keeping \(\sigma_R = 0.1\).
Show code
params_low_sigma_low <- params_low_sigma
params_low_sigma_low$h <- scenario_low$params$h
params_low_sigma_high <- params_low_sigma
params_low_sigma_high$h <- scenario_high$params$h
Fmsy_low_sigma <- calc_reference_points(params_low_sigma_low)$Fmsy
Fmsy_high_sigma <- calc_reference_points(params_low_sigma_high)$Fmsy
F_series_alt <- ifelse(block_id_low %% 2 == 1, Fmsy_high_sigma, Fmsy_low_sigma)
log_rps_list_fmsy <- vector("list", n_sims_alt_low)
for (i in 1:n_sims_alt_low) {
sim_alt_fmsy <- run_simulation_variable_h(
params_low_sigma,
h_series = h_series_low,
F_series = F_series_alt,
seed = 6000 + i
)
rps_fmsy <- sim_alt_fmsy$Recruitment[-1] / pmax(sim_alt_fmsy$SSB[-length(sim_alt_fmsy$SSB)], 1e-12)
log_rps_list_fmsy[[i]] <- data.frame(
Year = 2:n_years_alt_low,
Sim = i,
Log_R_SSB = log(rps_fmsy)
)
}
log_rps_fmsy_df <- do.call(rbind, log_rps_list_fmsy)
median_log_rps_fmsy <- aggregate(Log_R_SSB ~ Year, data = log_rps_fmsy_df, FUN = median)
p_alt_fmsy <- ggplot(log_rps_fmsy_df, aes(x = Year, y = Log_R_SSB, group = Sim)) +
geom_line(alpha = 0.2, size = 0.4, color = "steelblue") +
geom_line(
data = median_log_rps_fmsy,
aes(x = Year, y = Log_R_SSB),
inherit.aes = FALSE,
color = "black",
size = 1.1
) +
labs(
x = "Year",
y = "log(Recruitment / SSB)",
title = "Alternating steepness, regime-specific Fmsy (sigma_R = 0.1)"
) +
ggthemes::theme_few()
p_alt_fmsy
Even with a known alternation in steepness, detecting a productivity shift in practice can be difficult. Recruitment variability, observation error in SSB, and the lagged relationship between SSB and recruitment can mask regime changes, especially over short time windows. This makes it challenging to distinguish a true shift in productivity from stochastic swings, and underscores the value of trigger rules that rely on sustained signals rather than single-year anomalies.
9 Discussion
This comparison of low (\(h = 0.6\)) and high (\(h = 0.9\)) steepness scenarios demonstrates the critical influence of stock-recruitment steepness on fishery management. The main evidence comes from the reference points in Table 1 and the tradeoff patterns in Figure 1, Figure 2, and Figure 4. This simulation is a simplified, pollock-like evaluation of SPR-based control rules, grounded in the Tier 3 framework described in the BSAI FMP (North Pacific Fishery Management Council 2024). In this model, the proxy \(F_{MSY}\) is \(F_{40\%}\), with a sloping rule that holds \(F_{40\%}\) above 0.4 \(B_0\) and linearly ramps to zero at 0.05 \(B_0\). That structure aligns with SPR proxy guidance and the literature on SPR-based reference points (North Pacific Fishery Management Council 2024; Szuwalski and Punt 2025), while remaining intentionally stylized to match the scenarios tested here.
The two steepness scenarios (\(h = 0.6\) and \(h = 0.9\)) serve as a tractable proxy for productivity shifts discussed in climate and recruitment literature (Spencer et al. 2019; Punt et al. 2024; Szuwalski et al. 2023). We do not explicitly model environmental covariates or the ACLIM framework, but those references motivate the sensitivity analysis and the interpretation of performance differences across productivity regimes (Hollowed et al. 2020).
Finally, the catch-capped HCR variant (1.45 Mt) is a stylized analogue to ecosystem-cap effects that can dampen catch variability and conserve biomass during high abundance (Holsman et al. 2020; Ianelli et al. 2011). While the cap value here is illustrative and not a direct re-creation of the 2 Mt cap, the direction and intent are consistent with the cited findings, making the comparison appropriate for the simplified model context.
9.1 Management Triggers and Objectives
This analysis can be extended to explicitly test management objectives and triggers. Objectives define the desired outcomes and can be evaluated directly with the performance metrics already computed here (mean catch, catch variability, and risk of low biomass). Triggers then specify when the baseline capped HCR should tighten or relax specific elements of the rule (e.g., a lower F target, a smaller cap, or a steeper ramp).
Examples of objectives that map cleanly onto the current outputs include:
- Maintain low risk: keep \(P(SSB < B_{MSY})\) below a chosen threshold at terminal years (see Table 4 and Table 5).
- Maintain catch stability: limit catch CV over the last 10 years (see Figure 9).
- Maintain yield: keep mean catch above a target while satisfying the risk and stability bounds.
Illustrative triggers that could be implemented within the capped HCR framework include:
- Biomass trigger: if a 3-year moving average of \(SSB/B_{MSY}\) falls below 1.0, reduce \(F_{40\%}\) by a fixed fraction for the next 5 years.
- Trend trigger: if SSB declines for 4 of 5 consecutive years, lower the cap and reduce the ramp slope until SSB recovers above the threshold.
- Catch variability trigger: if catch CV exceeds a policy limit, temporarily tighten the cap or reduce the \(F_{40\%}\) target.
- Recruitment anomaly trigger: if recruitment falls below a historical percentile for multiple years, tighten the ramp parameters while maintaining the cap.
- Prolonged low recruitment trigger: if recruitment is below its long-term mean for 5 consecutive years, reduce the \(F_{40\%}\) target and lower the cap until recruitment rebounds.
These triggers can be evaluated by their false-alarm rate, time-to-recovery, and tradeoffs against yield, allowing the Council to align management actions with explicit objectives rather than relying on implicit assumptions about productivity.
Illustrative trigger table (capped HCR adjustments):
| Trigger | Metric | Threshold | Capped HCR adjustment |
|---|---|---|---|
| Biomass safeguard | 3-year mean \(SSB/B_{MSY}\) | < 1.0 | Reduce \(F_{40\%}\) target by 20% for 5 years |
| Persistent decline | SSB trend | Declines in 4 of 5 years | Lower cap by 15% and steepen the ramp |
| Catch instability | Catch CV (last 10 years) | > 0.30 | Temporarily tighten cap by 10% |
| Recruitment anomaly | Recruitment percentile | < 20th percentile for 3 of 5 years | Reduce \(F_{40\%}\) by 15% and tighten ramp |
| Prolonged low recruitment | Recruitment vs mean | Below mean for 5 consecutive years | Reduce \(F_{40\%}\) by 10% and lower cap by 10% |
9.2 Reference Points
| Metric | h = 0.6 | h = 0.9 |
|---|---|---|
| \(F_{MSY}\) | 0.352 | 1.28 |
| \(B_{MSY}\) | 2.66 Mt | 1.59 Mt |
| MSY | 1.31 Mt | 2 Mt |
| \(B_{MSY}/SSB_0\) | 0.32 | 0.19 |
9.3 Key Findings
Higher steepness dramatically increases sustainable fishing mortality: \(F_{MSY}\) increases from 0.35 to 1.28 (a 264% increase).
MSY is 53% higher with high steepness (2 vs 1.31 million tonnes).
\(B_{MSY}\) is lower with high steepness because recruitment remains high even at reduced biomass levels.
Management implications: Uncertainty in steepness translates directly to uncertainty in sustainable harvest levels. A stock managed as if \(h = 0.9\) when true \(h = 0.6\) would be at high risk of overfishing.
F40% sloping HCR vs constant F: The proxy Fmsy (F40%) with the sloping control rule yields a different catch/SSB tradeoff than the constant-F multipliers, and its mean realized F/Fmsy and outcomes shift between \(h = 0.6\) and \(h = 0.9\) (see Figure 7 and Figure 8).
10 References
Hollowed, A. B., J. N. Ianelli, P. A. Livingston, B. Planque, and W. S. Wooster. 2020. “Changed States of the North Pacific Ecosystem.” ICES Journal of Marine Science 77: 886–97.
Holsman, K. K., J. N. Ianelli, K. Aydin, and S. Zador. 2020. “Ecosystem-Based Fisheries Management Requires Ecosystem-Based Analysis.” Frontiers in Marine Science 7: 594.
Ianelli, J. N., K. K. Holsman, K. Aydin, and A. E. Punt. 2011. “Multi-Model Inference of Linkages Between Walleye Pollock and Pacific Cod Dynamics Across the North Pacific.” Marine and Coastal Fisheries 3 (1): 110–28.
North Pacific Fishery Management Council. 2024. Fishery Management Plan for Groundfish of the Bering Sea and Aleutian Islands Management Area. North Pacific Fishery Management Council.
Punt, A. E., C. S. Szuwalski, J. N. Ianelli, and A. B. Hollowed. 2024. “Evaluating Harvest Control Rules Under Ecosystem Change.” ICES Journal of Marine Science 81 (1): 45–62.
Spencer, P. D., J. N. Ianelli, A. B. Hollowed, A. E. Punt, and J. K. Parrish. 2019. “Modelling Spatiotemporal Dynamics of Walleye Pollock in the Bering Sea: A Long-Term Outlook.” Fisheries Oceanography 28 (2): 189–204.
Szuwalski, C. S., J. N. Ianelli, A. E. Punt, and P. D. Spencer. 2023. “Climate-Informed Recruitment Dynamics: Implications for Harvest Control Rules.” Fisheries Research 265: 106745.
Szuwalski, C. S., and A. E. Punt. 2025. “Steepness Uncertainty and Reference Points in Multistock Management.” Canadian Journal of Fisheries and Aquatic Sciences.
Appendix: R Session Info
Session information is reported below to document the computational environment used to generate the figures and tables.
Show code
sessionInfo()R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.2
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Los_Angeles
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggthemes_5.2.0 ggplot2_4.0.1
loaded via a namespace (and not attached):
[1] Matrix_1.7-4 gtable_0.3.6 jsonlite_2.0.0 dplyr_1.1.4
[5] compiler_4.5.2 tidyselect_1.2.1 stringr_1.6.0 splines_4.5.2
[9] scales_1.4.0 yaml_2.3.12 fastmap_1.2.0 lattice_0.22-7
[13] R6_2.6.1 labeling_0.4.3 generics_0.1.4 knitr_1.50
[17] htmlwidgets_1.6.4 tibble_3.3.0 pillar_1.11.1 RColorBrewer_1.1-3
[21] rlang_1.1.6 stringi_1.8.7 xfun_0.55 S7_0.2.1
[25] cli_3.6.5 withr_3.0.2 magrittr_2.0.4 mgcv_1.9-3
[29] digest_0.6.39 grid_4.5.2 nlme_3.1-168 lifecycle_1.0.4
[33] vctrs_0.6.5 evaluate_1.0.5 glue_1.8.0 farver_2.1.2
[37] codetools_0.2-20 rmarkdown_2.30 purrr_1.2.0 tools_4.5.2
[41] pkgconfig_2.0.3 htmltools_0.5.9